{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating e-ph coupling matrix element (VASP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we calculate the overlap `<ψ_i(0)|ψ_f(Q)>` and its derivative.  \n",
    "`ψ_i(0)` is a perturbed band-edge state, and `ψ_f(Q)` is a localized defect state at the configuration coordinate `Q` which can be calculated using `get_del_Q.py`.\n",
    "The electron-phonon coupling matrix element is given by   \n",
    "`W_if = (ε_f - ε_i) d <ψ_i(0)|ψ_f(Q)> / dQ`.\n",
    "\n",
    "For details, see Eq.13 in Ref 1. \n",
    "\n",
    "We use [`pymatgen`](https://kylebystrom.github.io/pawpyseed/) and [`pawpyseed`](https://pymatgen.org) to calculate the overlap of **all** electron wave functions.\n",
    "\n",
    "[1]\tA. Alkauskas, Q. Yan, and C. G. Van de Walle, Phys. Rev. B 90, 075202 (2014)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overlap\n",
    "\n",
    "Use `get_wf.py` in `/script` directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "$ python get_wf.py   -d <band_index of ψ_i> -b <band_index of ψ_f>  -D <path to the initial geometry> -i <path to the deformed geometry 1> [<path to the deformed geometry 2> ...]\n",
    "\n",
    "GRID ENCUT 918.2911873931497\n",
    "finished making projector list\n",
    "--------------\n",
    "ran get_projector_list in 0.033210 seconds\n",
    "---------------\n",
    "STARTING PROJSETUP\n",
    "started setup_proj\n",
    "calculating projector_values\n",
    "onto_projector calcs\n",
    "Done\n",
    "--------------\n",
    "\n",
    "...\n",
    "\n",
    "=================================\n",
    "----------- Overlaps ------------\n",
    "=================================\n",
    "[[0.0574]\n",
    " [0.0267]\n",
    " [0.0345]\n",
    " [0.0643]]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## W_if calculation\n",
    "\n",
    "Use `get_del_Q.py` to get a distance between structures (`ΔQ`) in `amu^0.5Å`. `ΔQ` is a collective variable of the mass-weighted deformation `ΔQ = ΣmΔ𝐑`.\n",
    "\n",
    "`ϵ_i` and `ϵ_f` are KS eigenvalues of the initial (`ψ_i(0)`) and final states (`ψ_f(0)`).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "W_if_n = 0.13949890444622062\n"
     ]
    }
   ],
   "source": [
    "# using Plots\n",
    "# \n",
    "ΔQ = 0.18059526482789678 # ΔQ between DISP_000 and DISP_001 in [amu^0.5Å]\n",
    "ϵ_f = 5.405963  # unoccupied trap\n",
    "ϵ_i = 6.227692  # CBM\n",
    "\n",
    "overlap = [0.0574, 0.0267, 0, 0.0345, 0.0643]\n",
    "overlap[1:3] = -overlap[1:3]\n",
    "overlap[3] = 0\n",
    "# # Central finite difference\n",
    "diff = [1/12 -2/3 0 2/3 -1/12] * overlap / ΔQ\n",
    "\n",
    "W_if_n = (ϵ_i-ϵ_f) * diff[1]\n",
    " \n",
    "println(\"W_if_n = \", W_if_n)\n",
    "# plot(overlap)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.1",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
